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Abstract 

A perturbation theory of the static response of insulating crystals to homoge- 
neous electric fields, that combines the modern theory of polarization with the 
variation-perturbation framework is developed, at unrestricted order of per- 
turbation. Conceptual issues involved in the definition of such a perturbative 
approach are addressed. In particular, we argue that the position operator, 
X, can be substituted by the derivative with respect to the wavevector, d/dk, 
in the definition of an electric-field-dependent energy functional for periodic 
systems. Moreover, due to the unbound nature of the perturbation, a regular- 
ization of the Berry-phase expression for the polarization is needed in order to 
define a numerically-stable variational procedure. Regularization is achieved 
by means of a discretization procedure, which can be performed either before 
or after the perturbation expansion. We compare the two possibilities, show 
that they are both valid, and analyze their behavior when applied to a model 
tight-binding Hamiltonian. Lowest-order as well as generic formulas are pre- 
sented for the derivatives of the total energy, the normalization condition, the 



eigenequation, and the Lagrange parameters. 
77.22.-d,77.22.Ch,78.20.Bh,42.70.Mp 
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I. INTRODUCTION 



Until the early nineties, the formulation of a quantum-mechanical theoretical framework 
for the study of the physics of electric polarization in solids had remained a challenging 
problem.011 Even the definition of the polarization itself as a bulk quantity, independent of 
surface termination, was the subject of heated debate.i^ 

This picture changed when King-Smith and Vanderbilt (KS-V)il proposed a formulation 
(the modern theory of polarization - MTP), which resolved the conceptual difficulties as- 
sociated with the definition of this quantity for continuous, periodic, charge distributions. 
In their work, the electric polarization of an insulating crystal is related to a Berry phaseB 
computed from the valence wavefunctions. The existence of a band-structure Berry phase 
had already been discussed by Zak and coworkersj^ before its connection with the electronic 
polarization was established by KS-V. Besides settling the important conceptual question 
related to the definition of polarization as a bulk quantity, the KS-V theory provided an 
entirely new framework for the computation of the polarization of a crystal maintained at 
vanishing homogeneous electric field. Since its formulation, the theory has been examined 
in greater detail by KS- V0 and by Rest J, and extended to many-body systems by Ortiz 
and MartinO. The relation between polarization and the phases of the wavefunctions has 
also led to a reexamination of the role this quantity plays in the Density Functional Theory 
(DFT)ii formulation of the ground-state properties of extended systems.li3ll3 

Of no less importance is the conceptual relationship between the spontaneous polariza- 
tion and the centers of charge of the Wannier functions (WF) of the occupied bands, which 
was also discussed by KS-Vi and previously by Zak.H This connection was later generalized 
by Nunes and Vanderbilti^ (NV) to deal with insulators placed in non-zero external homo- 
geneous electric fields: they introduced field-dependent "polarized" WF's and a method for 
their computation. NV argued that, in the static-response regime, the state of an insulator 
under an external homogeneous electric field is one in which the periodicity of the charge 
density is retained, despite the fact that the perturbation lacks the lattice-translational 
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symmetry of the unperturbed crystal. Such state is actually a long-lived resonance of the 
system, as rigorously demonstrated by Nenciul^. 

It is well known that within DFT0lii, ground-state properties of condensed-matter sys- 
tems, such as equilibrium lattice parameters, bond lengths and bond angles, among others, 
can be obtained with an accuracy of a few percent in comparison with experimental results. 
Within DFlQIil, the NV method has recently been applied to the computation of the polar- 
ized WF's and the dielectric constant of silicon and gallium arsenideil. The latter quantity 
is related to the change of polarization due to a change of homogeneous electric field, in the 
linear regime, or equivalently, to the second-derivative of the total energy with respect to 
the homogeneous electric field. 

Specific treatments have been developed (noticeably within DFT) for the study of the 
response of crystals to "external" perturbations, like phonons, stresses or homogeneous 
electric fields. The latter, on which we focus exclusively, can be taken as a homogeneous 
field or as the limit of long-wavelength perturbations. 

In the so-called direct approach!!, supercell calculations are employed to study both 
the unperturbed and perturbed systems, with the response functions being obtained by nu- 
merical finite-difference analysis of the changes induced by a long-wavelength perturbation 
applied to elongated supercells. The non-linear response regime is directly accessible, al- 
though it must be disentangled from the linear response of the system. However, because 
of the use of supercells, the computational cost of this approach is rather high. 

Alternatively, the specific response to a homogeneous electric field was considered within 
perturbation theory, already in the sixties. In the Random-Phase Approximation^ (no local- 
field effectH included), the response of the wavefunctions is obtained through a sum over 
states, involving matrix elements of the position operator between valence and conduction 
stateJ^. This technique was generalized to the computation of second- and third-order 
susceptibilities^. The need to compute many unoccupied bands is the bottleneck of this 
method. 

Local-field effects can be reintroduced on top of such a sum-over-states approach either 



in a matrix-inversion frameworkilii, or in an iterative approach!^. In their calculation of 
linear susceptibilities, Levine and Allan included a so-called "scissor-operator" correction, 
that was understood later to compensate some deficiencies of local-density approximation 
computation of long-wavelength response functionJi^. Also, Levine and coworkers proposed 
rather involved expressions for the second-order and third-order DFT susceptibilitiesi^. 

To a large extent. Density Functional Perturbation Theory (DFPT)&li overcome the 
limitations of the previously mentioned approaches, at the price of non-negligible additional 
coding. At the lowest order in the homogeneous-electric-field perturbation, this method was 
introduced by Baroni and collaboratorsil'ii. It is based on an iterative solution for the first- 
order change in the wavefunctions, which allows for the self-consistent inclusion of local-field 
effects, besides eliminating the cumbersome sum over conduction bands. It does not employ 
supercells, and can be applied to perturbations of arbitrary wavelengths. In the DFPT, the 
computational workload involved in the computation of linear-response functions is of the 
same order of that involved in one ground-state calculation. 

DFPT is part of a class of formalisms in which perturbation theory is applied to a 
variational principleS. This interesting combination leads to a generic "2n+l" theoremil0, 
as well as variational properties of even-order derivatives of the energ For example, 

one can compute the third-order derivative of the energy from the first-order derivative of 
wavefunctions, and the fourth-order derivative of the energy is variational with respect to 
the second-order derivative of wavefunctions. The expressions derived in this framework are 
surprisingly simple and can be formulated at all order of perturbations. 

However, the treatment of homogeneous electric fields in this variation-perturbation 
framework is plagued by difficulties similar to those encountered in the theory of polar- 
ization. Shortly after the appearance of the MTP, Dal Corso and MauriS, building upon 
the NV work, proposed a very concise expression for the second-order susceptibility which 
was later applied successfully to compute this quantity for a variety of systems^. 

In the present work, we formulate a perturbation theory of the static response of insulat- 
ing crystals to homogeneous electric fields that combines the conceptual ideas of the MTP 



with the variation-perturbation theory. A major achievement of our work is the presenta- 
tion of formulas vahd for unrestricted order of perturbation theory. We also examine the 
low-order expressions in some detail, and recover the expression proposed by Dal Corso and 
Mauri0. 

The theory is worked out directly in reciprocal space, in terms of the Berry phase as- 
sociated with the occupied bands of the perturbed crystal, in the manner of the MTP. 
The conceptual issues involved in the definition of a perturbative approach for the problem 
are addressed. The Berry phase is argued to remain a valid concept in the presence of the 
periodicity-breaking electric field. The periodicity of the charge density is assumed to survive 
the application of the field, and the Berry phase is obtained from the associated polarized pe- 
riodic wavefunctions. By working out the perturbative approach in terms of these polarized 
states, we obtain very compact expressions for the high-order dielectric-response functions 
of the crystal. These can be numerically obtained on the basis of iterative equations for the 
second- and higher-order terms in the perturbation expansion of the wavefunctions of the 
system, as in the DFPT approach for other types of perturbation. We will not deal explicitly 
with the exchange and correlation parts of the DFT functional: the main difficulties that we 
want to address in this paper are not related to them. The formalism can be extended to 
include exchange and correlation terms in a self-consistent fashion, in the manner presented 
in Refs 0,^. 

Any application of the MTP involves a discretization of the Berry-phase expression, in 
terms of a series of wavevectors for the electronic wavefunctions. We have discovered that 
such a discretization, that appears naturally also in the present framework, can be performed 
at two different conceptual levels when merged with perturbation theory: either after the 
derivation of formal expressions at different orders of perturbation theory, starting from a 
continuous Hamiltonian, or at the level of the field-dependent Hamiltonian itself, before any 
perturbation expansion is performed. We will refer to the first approach as the discretization 
after perturbation expansion (DAPE) formulation, and to the second as the perturbation 
expansion after discretization (PEAD) formulation. 
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In order to judge the relative merits of these two approaches, and also the correctness 
of the global framework, we analyze the behavior of a model tight-binding Hamiltonian, for 
which analytical responses to an electric-field perturbation have been obtained up to the 
fourth order. 

The paper is organized as follows. In the next section, we address the conceptual is- 
sues associated with the definition of a Hamiltonian and its perturbative expansion for the 



electric-field problem. Sec. |T| summarizes the main results of the variational perturbation 



theory which are used in this work. In Sec. we work out the continuous formulation of 
the problem and its perturbation expansion, from which we obtain the DAPE version of our 
theory. In Sec. 0, we work from the start using a discretized expression for the polariza- 
tion, which leads to the alternative PEAD formulation. The theory is applied to a model 



one-dimensional system in Sec. VI. 



II. INSULATORS IN AN ELECTRIC FIELD: CONCEPTUAL 

CONSIDERATIONS 

A. The modern theory of polarization 

In the MTpi'0il, the change in electric polarization per unit volume induced by an 
adiabatic change in the crystalline potential (the self-consistent Kohn-Sham potential in the 
context of DFT) is written 

AP = £^dX = P{X,)-P{X,), (1) 

with P(A) given in terms of a Berry phase associated with the occupied bands of the system 

«S) , (2) 

where — e is the electron charge, A is a parameter representing the adiabatic change in the 
potential, and the factor of 2 in the numerator accounts for spin (in this work we will consider 
only spin-unpolarized systems). 
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The gauge relation between periodic functions Mnk+G(r) = e~'*^''"Mrik(r) is established by 
requiring that the Bloch eigenstates be periodic in reciprocal space, i.e. ipnk = V'nk+G; where 
G is a reciprocal-lattice vector. With this choice of gauge, the polarization changes given 
by Eq. |I] are defined to within a factor of (2e/f2)R, where R is a lattice vector. i Eq. ^ was 
derived under the restriction that the macroscopic electric field inside the crystal vanishes. 
Moreover, it also requires that the set of wavefunctions be different iable with respect to k.@ 

The actual evaluation of the polarization in Eq. |^ is carried out on a discrete mesh of 
points in reciprocal space. Because this expression depends on the phase relationships of 
wavefunctions at different k-points, the following discretized version was proposed by KS-V: 

^11 = y {in det [ ( u% I u%^^ )] } ; (3) 

where P\\ is the component of the polarization along the direction of a short reciprocal- 
lattice vector, G||, and A^^ is the number of k-points sampling the Brillouin zone along that 
direction for each value of k_|_, with kj = k^ + jG\\/Nk. 

From a calculational point of view, this discretized expression ensures that the final result 
is unaffected by random numerical phases which may be introduced in the wave functions at 
different k-points, when these are independently determined by the diagonalization routine. 
However, Resta has taken the view that the discretized expression is to be regarded as more 
fundamental than the continuous form.0 For the formulation of the electric-field response 
that we develop in the present work, discretization is crucial in order to define a numerically- 
stable minimization procedure. We will come back to this point in Sec. |TD. 



The Berry-phase expression can be transformed into a real-space integral involving the 
Wannier functions of the occupied bands, leading to a physically-transparent expression for 
the polarization in terms of the centers-of-charge of the Wannier functions:i'iHl§'@ 

P(A) = -|E/rhi'^f^r, (4) 

where Q is the unit-cell volume. 

In principle, the above expressions are valid only at vanishing electric field. However, 
it was soon reahzedSill that Eq. ^ could be extended to the non-zero field problem, by 
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introducing the so-called polarized Wannier functions. Polarization effects were then related 
to the field-induced shifts in the Wannier-f unction centers of charge. In the present work, 
Eq. is argued to apply to the non-zero-field problem, thus defining a field-dependent Berry 
phase containing the polarization effects. This allows us to work out a perturbative approach 
for the finite-field problem. 



B. Definition of the energy 

The study of the problem of insulators under an external electric field has traditionally 
met with conceptual difficulties, related to the non-analyticity of the perturbation, as dis- 
cussed in detail by Nenciu0. Upon the application of the external field, the spectrum of 
electronic states changes non-analytically, with the band structure of the insulator at zero 
field being replaced by a continuum of eigenvalues spanning the entire energy axis (from 
— oo to +00), even for a field of infinitesimal strength. From a mathematical point of view, 
the unbound nature of the perturbation term, eS-r (hereafter, we use S for the magnitude of 
the electric field), hinders the straightforward application of perturbation theory, since the 
diagonal elements of the position operator in the basis of the unperturbed Bloch states are 
ill-defined§§'0. Strictly speaking, an infinite crystal in the presence of an external electric 
field does not have a ground stateElS. 

From the physical point of view, in the limit of weak to moderate fields, the tunneling 
currents (which destroy the insulating state at sufficiently high fields) can be neglected, and 
only the polarization of the electronic states by the external field is considered. One is left 
with a picture of the problem in which the insulating state of the unperturbed crystal is 
preserved, hence the band structure and the periodicity of the charge density are retained. 
The problem is thus physically well defined. The theory we develop below will concern this 
periodic polarized insulating state, resulting from the application of the electric field. Avron 
and Zak0 have discussed the stability of the band structure in the present context, while, as 
mentioned in the introduction, Nenciu0 has rigorously shown such state to be a long-lived 
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resonance of the problem, in the regime of weak fields. 

Building on the conceptual framework set by the MTP, NV introduced a practical real- 
space method to handle the problem on the basis of the WF formulation of the polarization. 
They noted that this polarized insulating state can be represented in real space by a set 
of field-dependent polarized WF's. In this way, the relationship between the WF's centers 
of charge and polarization can be extended to the non-zero field situation, and an energy 
functional is defined as follows 

E [{w^}, S] = [{w^]\ - n S P [{w^}] , (5) 

where {w^} is the set of field-dependent Wannier functions. 

Because the state underlying the above expression is not a ground-state, rather a reso- 
nance, the energy functional is only well defined for WF's of finite range.0 The truncation 
of the WF's provides a mathematical procedure for the regularization of the problem. We 
will come back to this pathology in section [11 L)| 

As a corollary of the existence of polarized WF's, we consider now the representation of 
the system in terms of polarized Bloch orbitals. In the following sections, we develop two al- 
ternative formulations in which the A;-space MTP expressions for the electronic polarization, 
Eqs. I and |, are extended to the non- zero-field problem. We will also have to regularize 
the ensuing expressions, which is possible by means of discretization of the /c-space integrals 
in such a way that polarized valence bands are stable against mixing with the conduction 
bands. 

For simplicity, we concentrate on a one-dimensional non-interacting spin-unpolarized 
system. The Hamiltonian for the unperturbed periodic insulator is given by 

if(o) =k+Vq, (6) 

where K is the kinetic-energy operator, and Vq is a periodic potential, i.e., [Vq, Ti] = 0, where 
Ti denotes a translation by a lattice vector. The zero-field function u^^l is the periodic part 



of the unperturbed Bloch-orbital, obeying the eigenvalue equation Hj 
where 
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is the unperturbed cell-periodic Hamiltonian. 

Under the action of an external electric field, the Hamiltonian becomes 

H = + eSx . 



(8) 



As discussed above, this Hamiltonian is not amenable to a conventional perturbation treat- 
ment. We can arrive at an expression that is applicable to extended periodic systems, from 
the following considerations. 

First, let us assume that we can define a set of field-dependent cell-periodic functions, 
representing the polarized state of the system. These are the Fourier transform of the field- 
dependent Wannier functions introduced by NV. Eq. |^ is thus extended to the non-zero field 
problem, with 



nk 



d_ 
dk 



(9) 



defining a field-dependent polarization, including the spontaneous and induced parts of this 
quantity. 

Next, combining this definition with Eq. |^ we write 
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u 



nk 



(0) 



u 



nk 



+ s(u: 



It 



dk 



u 



nk 



(10) 



for the total energy in terms of field-dependent cell-periodic functions, with the unit-cell 
volume VL = a for our ID system. 

Consider now the expansion of the set {u^k\ terms of the complete set of zero-field 
periodic functions: 



11^ 
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In terms of this expansion, Eq. |T0| is written 

N oo 
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■ (12) 
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This expression is in the form of the matrix representation, in the basis {timi}, of the 
expectation value of the following "operator" : 



This suggests Eq. as an ansatz for the cell-periodic Hamiltonian, now including the 
perturbation term f/|"^'^* = ieS-^. 

In subsection [11 Q we show that the perturbation Hamiltonian operator given in Eq. ^ 
by its matrix representation, 



is the periodic part of the x operator. 

Another way of arriving at this ansatz is by the following argument. Consider the first- 
order change in the total energy which can be obtained, without postulating the existence 
of the field-dependent functions, by combining Eqs. | and ||, as follows: 



This is simply the coupling of the spontaneous polarization of the system to the external 
field. From textbook perturbation theory, the first order change in the energy is given by 
the diagonal matrix elements of the perturbation term, leading again to the form of the 
Hamiltonian in Eq. |T3|. 

Some remarks are needed about the application of perturbation theory to this Hamil- 
tonian. Strictly speaking, Umn{k) is not an operator (its transformation properties under 
unitary transformations will be discussed in subsection |11 C] ). However, it can be shownEl 
that this is reflected only in the first-order change of the single-particle eigenvalues, when 
we analyze single-particle quantities obtained in the perturbation expansion of Eq. |T3|. All 
other single-particle quantities, such as wave-function derivatives and higher-order eigen- 
value derivatives are actually gauge- invariant, and thus well defined. More importantly, in 
the following developments we will be interested only in quantities that are integrated over 




(13) 




(14) 




(15) 
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the Brillouin zone (the derivatives of the total-energy with respect to the apphed field), 
which will be shown to be gauge invariant. As a final observation, it can also be shownS 
that under unitary transformations, the first-order eigenvalue acquires a change which, when 
integrated over the Brillouin zone, leads to a first-order energy derivative which is defined 
modulo the quantity —eSi, where i = Na is a lattice vector {N is an integer). This is 
consistent with the fact that, in the MTP, the zero-field polarization itself is defined modulo 
-ei. 

The continuous formulation of our theory is based on the application of a variational 
perturbation treatment to Eqs. [l^ and [13|. Alternatively, the PEAD formulation is derived 
by applying the variational principle to the total energy written in terms of the discretized 
form of the polarization. In this case, we combine Eqs. ^ and |^ to write 



E 



{Unk,};£] = ^ ] ^ ^(Wnfc, |i^i°^|Mnfc,) - ^I]lm{lndet [5™„(/Cj-,/Cj+l)]} I , (16) 



N TVfe 



e£ 



n=lj=l j=l 



where j runs over the Nk /c- vectors in the discretized Brillouin zone, Ak = 271 /aNk, and 

(17) 

is the overlap matrix between states at adjacent points in the reciprocal-space mesh. 



C. Position operator for periodic systems 

In view of the above discussion, we examine now the action of the position operator in 
a space of periodic functions. Keeping in mind that we wish to retain the periodicity of 
the charge density, we seek to arrive at a consistent definition for the action of x in that 
space. This problem has been recently tackled by RestaEl, who suggested an intrinsically 
many-body redefinition of x, in the context of periodic systems. In the spirit of retaining 
a single-particle picture, here we only offer a heuristic justification for the form of the 
perturbation term given in Eq. [T^. For this, we use the crystal momentum representation 
(CMR), following the discussion in the paper by Blount.@ 
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Let f{x) denote a square- integrable function. The full set of zero-field Bloch eigenstates 
of a periodic Hamiltonian forms a complete basis to expand f{x): 

m = ^J dkY^^^2{x)Uk) = ^J dke^'^Y.^^^lix)Uk) . (18) 



The action of x on f{x) is given in the CMR by the expression 

.dfnik) 



xf{x) = ^ j dke^^^Y.^f,{x) 
Ztc J „ 



+ J2^nn'ik)Uik) 



(19) 



where Unn' is defined in Eq. 0. 

Blount examined the transformation properties of the two terms appearing in Eq. |l^, 
with respect to the choice of the phases of the Bloch orbitals. Consider the following CMR 
decomposition of x: 

d 

X = i—+Unn'ik) = Xd+Unn'{k) , (20) 

where Xd = i-§j^ is diagonal in the band index. He showed that when the Bloch orbitals are 
multiplied by a phase factor e*'^"*^'^\ the term Xd transforms as x'd = Xd — 5nn'd(pnik) /dk, 
while a compensatory change occurs in the diagonal term Unn- So, the two terms in Eq. ^ 
do not transform separately like operators, while their sum does. In our formulation, we 
use the second term on the right, Unn', to define a periodic Hamiltonian for the electric-field 
problem which, from this discussion, is not by itself an operator in the strict sense. 

We show now that Unn'{k) is translationally invariant, while Xd (like x) is not. Consider 
a translation Ti by a lattice vector ^. The commutation relation for x and Ti is written. 

[x,T,]=£T,. (21) 

To obtain the commutation relation of the perturbation term in Eq. 0, we expand 
Eq. (pl|) in the CMR representation. Let g{x) = Tif{x) = f{x — i). From Eq. (p!8| ) and 

Unkix - ^) = wife (a;) we get 

gix) = dke^'^^^-'^j:^^^lix)Uk) . (22) 
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From this expression, it follows that 



9n{k) = e-'''Uk). 



(23) 



Further, from Eq. ( |T9|) we obtain 
while from Eqs. (O) and (ESf) we get 



. dfnjk) 
dk 



+ Y.Unn'{k)U,{k) 



(24) 



xT,f{x) = xgix) = ^ / dke''^J2-^'^ 

= ^jdke^'^^-'^Y.^^l{x 



.dgnjk) 
' dk 



+ ^Unn'{k)gn'{k) 



ifnik) + + Y.^nn'{k)fAk) 

ok , 



(25) 



where, in the last step, we use the result i-^gn{k) = e'"^^^ Unik) +i-^fn{k) . Combining 
Eqs. (24) and (25) we arrive at the CMR expansion of the commutation relation in Eq. (^TJ). 
Moreover, the above development immediately shows that 



[Unnik),Te]=0. 



(26) 



From the above, we observe that the perturbation term in Eq. [T^ is invariant under 
lattice translations. 



D. Wannier-function cutoff in real space 

The definition of an energy functional at finite fields requires careful analysis. Because 
the problem does not have a ground state, a regularization procedure is required for the 
definition of a numerically-stable functional capturing the physics of the state of the system 
after the electric field is turned on. 

In the NV treatment of the problem in real-space, regularization is achieved with the in- 
troduction of truncated Wannier functions, which are constrained to zero beyond a real-space 
cutoff -RcS As discussed by NV, in the limit ^ oo the functional becomes pathological, 
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with the property that a state having an arbitrary value for the polarization can be con- 
structed without changing the value of the energy, when working at fixed polarization, or 
conversely with the development of a growing (infinitely many in the Rc oo limit) false 
local minima when working at fixed electric field. 

In order to develop this analysis on a sound mathematical basis, one performs a Legen- 
dre transformations, from the £^-dependent total energy E[£] to the P-dependent electric 
enthalpy E[P]: 

E[P] = inf {E[S] + aPS} . (27) 
The total energy was obtained previously [see Eq.(5)] thanks to the trial Wannier functions, 
E[S] = inf {E[{w}; S]} = inf {e(^\{w}] - aSP[{w}]} , (28) 
while a constrained search alternatively gives its Legendre transform, 

|?«}sucnthatF[{ui|J=_P ^ 

The zero-electric field total energy functional of the Wannier functions is -E'*^'^^[{w}] = 



We aim at understanding the pathologies of E[S] by examining its expression as the 
inverse Legendre transform of E [P] , 

E[S] = inf [e[P] - aPS} , (30) 

for which we need to characterize the minima of -E[-P], as well as their local behavior. 

We consider, for simplicity, the case of a single occupied band. For a given finite value 
of Rc, the electric enthalpy is a periodic function of P, and E[Pq] = E^^^[{wo}] = Eq is the 
zero-field ground-state energy {wq is the zero-field valence-band Wannier function). For large 
values of Rc, it becomes possible to build a set of ^-dependent functions (to be normalized), 

\w)= \wo) + p'/'r'/' \w^\e)) , (31) 
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with arbitrary value of P, where wq is a zero-field valence-band Wannier function centered 
at the origin, and WQ^{i) is an empty conduction-band function centered at the site i within 
the range of Rc, whose coefficient is on the order of We consider the lattice constant 

a to be the unit of length. The energy for these states is 

Eo + Pr'{Ecb-Eo) , (32) 

where Feb is the expectation value of the energy for the conduction band Wannier function. 
Due to the exponential decay of Wannier functions for insulators, the value of the polarization 
is 



(w 








Wo) + pr^ {wf{i) 


i7(o) 
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w \w 


) 




1 + P£-1 







{w 
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\w) 


{w 


w) 



Po + P, 



(33) 



since {w'^{t} 



X 



1 + P£-i 

'^o'i^)) — -\- Pq'- In the hmit £ — > oo, these wavefunctions have an 
arbitrary value of P, and an energy infinitesimally close to Eq. The E[P] curve becomes fiat 
in this limit, and only derivatives in an infinitesimal region around the ground-state solution 
remain well defined. The development of multiple minima at finite fields corresponds to the 
same situation, as a growing number of minima with energies that become degenerate in the 
Rc ^ oo are associated to states with different values of polarization. No global minimum 
as a function of polarization can be found. In the next subsection, we analyze the behavior 
of the energy functional for a model system in reciprocal space. We will show that the same 
pathology manifests itself in the limit Ak — > 0, where Ak is the discretization of the mesh 
of /c-points in the Brillouin zone. 



E. Reciprocal- space analysis of a model system 

For the present analysis, as well as for the application of the perturbation expansions to 
be developed in the following sections, we chose a one-dimensional (ID) two-site periodic 
model defined by two parameters, the hoping integral t, and the on-site term which we 
choose as —A/2 and A/2, for sites 1 and 2, respectively. The Hamiltonian can be rescaled 
by A to become a one-parameter ^ t) model, defined as 
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H = [^4,1^2,1 - ^c\ iCi^i + t [4 ,C2,/ + 4iCi,,+i + /i.e.] I , (34) 

where I runs over unit cells. Whenever we are concerned with the ID model, we will consider 
all distances to be rescaled by the unit-cell period (i.e. we set a = 1 in the present section, 
in Sec. pl[ and in Appendix 0), such that on each cell, denoted by the integer /, we have 
the basis functions 0i(/) and 02(^ + 1/2)- 

We apply Bloch's theorem to write the Schrodinger equation for the cell-periodic func- 
tions: 



(0) 



„.(0)\ _ (0) 



U 



(0) 

nk I 1 



(35) 



where is the zero-field cell-periodic Hamiltonian. In the basis of periodic functions 
Xi = E« 01 (0 and X2 = E/ 02(^ + 1/2) we have 



(0) 



^ -- 2tcos|^ 
2 ^ 



2t cos I 



The corresponding secular equation, det 
ues 



(36) 



0, is easily solved for the eigenval- 



-(0) 



\ + At^ cos^ ^ 
4 2 



k' 

1 + cos^ - 



(37) 



where A = At. Negative and positive eigenvalues correspond to valence and conduction 
bands, respectively. Because the Hamiltonian is real, we can use the following parameteri- 
zation for the corresponding eigenstates 



cos Qh ^ 
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where a^k and ack are real numbers, with no lack of generality. Coming back to the eigenvalue 



equation H), 
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u. 



-(0) 
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u^^l ), we obtain 
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2t cos I 
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Integrating Eq. |3^ over the Brillouin zone gives the energy per unit ceU: 



1 r'^'^ 



dke 
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vk 



2 rf 



dy 1 + cos^ y 



(39) 



(40) 
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In order to discuss the pathology of the finite electric-field functional in A;-space, we 

consider a set of trial cell-periodic functions 

/ 



\Uk) 



(41) 



sin 9fc e*'^*= j 

for k G [— 7r,7r], where 9^, a^, and Pk are real numbers. Imposing the condition \uk+G) = 
e**^^ ,i we obtain 

cos ek+2n 6^"'=+^^ = cos Gfc 6*"'= 

sin 6^+2. e*'''=+2- = - sin 6^ e'"'' , (42) 

which leads to — ao = NaTr and j32n — Po = Npn. 

The expectation value of the zero-field Hamiltonian in the set of trial wave-functions 
gives 

^^°^[K}] = - r (^k Uk) = - r dk -Jcos(20fc) + 2tcos ( | ) sin (26^) cos7fc 

TTJO ^ 'TTJO I 2, \Z J 

where Jk = ak-Pk- 

The polarization for the trial state is 



(43) 
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+ / dk cos (29^ 

JO 



dk 
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Minimization of -E'^^-' [{-Ufc}] with respect to 6^ and 7^, by setting dE^^ /dQk = and 



a4°V^7fc = 0, with E^' 



(0) 



Uk 



(0) 



Ukj, leads to 



tan(29fc) = — 4tcos ( — 1 cos7fc 



sin 7fc = 
or 

2tcos(|) sin (26^) = . 
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(45) 



(46) 



At cos(A;/2) = 0, the solution of Eqs. ^ and ^ leads to sin(20fc) = and also implies that 
7fc is undefined. At cos(A;/2) 7^ 0, a minimum solution is obtained by setting 

sin = ^ 7fc = A^^ X 27r , 

and (47) 
tan(26fc) = -4tcos (^| 

The ground-state solution is given by •jk = (i.e., = Pk as in Eq. |3^) for all values 
of k, with 9fc defined by Eq. Note that a solution where 7^ jumps by a multiple of 
27r at A; = j^Tr is also consistent with Eqs. ^5| - ^ , but not with the restriction that be 



different iable with respect to k. Note also that, due to inversion symmetry, the zero-field 
ground-state polarization must vanish (modulo — e). This is what is obtained from Eq. 
by setting d'jk/dk = 0, = Pk, and — ao = NaTc. 

We consider now a trial wavefunction where 0^ is the same as in the ground-state 
solution, while 7^ behaves as shown in Fig. where it jumps by a value of 2tt over an 
interval Ak centered at an arbitrary value of k. We show now that in the Ak limit 
this function can be tailored to give an arbitrary value of the polarization, while its energy 
differs from the ground-state by an infinitesimal amount, of order Ak. 

The change in polarization for this state, with respect to the ground-state solution, is 
given by 

AP = ^ J^^ dk cos (29^) ~ -e cos {201^^)) • (48) 

AP in the above equation assumes values between — e and — e/(l + 16t^)^/^. By adding 

another kink in the definition of 7^, where this function changes by — 27r, we can build a 

solution having any arbitrary value of P in the interval — e [0, 1]. 

Let us consider the change in energy of the trial state. The function cos7fc, as shown in 

Fig. ^, differs from one over a small interval of the order of Ak. The change in energy with 

respect to the ground state is then 

1 / k\ 2t f (k)\ / \ 
AE = -J dA;2tcos ( - j sin(20fc) (cos7fc - 1) ^ cos I ^ j sin (^2e(fc) j x AA;. (49) 
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So, for the trial state AE when Ak 0. The E^^\P) curve becomes flat in 
this limit. This is the same pathology as the one discussed by NV in the real-space case. 
Discretization of the A;-space mesh in the Brillouin zone is thus essential for the numerical 
stability of the energy functional. 

Now, we show that for the discretized version of the formulation, a change in P implies 
a finite change in the energy. The discretized polarization is written 



P 



^ Im In (uk 
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Uk 



'j+i 



- <^ ANp + Im ^ In 



TT 



k=l 



COS Qkj COS 6fc^^-^e*^'^''j+i '""j ' 



+ sin 6fc, sin 6/ 



'j+i 



(50) 



with the energy given by 



Uk, 



Nk 



H, 
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Uk 
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-\ COS (26fc^.) + 2tcos f y j sin (26^^.) cos7fcj. 



Again we consider the ground-state solution -^kj = 0, with 6^^. given by Eq. An 
arbitrary change in polarization can be introduced by setting 7^^. 7^ at a given kj, while 
keeping the values of G and 7 at all the other fc-points unchanged. In this case, it can 
be immediately seen that a finite change AEq = (At/Nk) cos (^y) sin (20A:j) (cos7a,.^. — 1) is 
introduced in the discretized energy. 



F. Summary 

The theoretical treatment of a periodic insulator placed in an homogeneous electric field 
is plagued by severe conceptual difficulties: (1) the potential associated with an electric field 
is non-periodic and unbounded; (2) for that reason the spectrum of electronic states changes 
non-analytically upon the application of a homogeneous electric field; (3) the quantity con- 
jugated to the electric field, namely the polarization, cannot be computed as the expectation 
value of the position (or any other) operator; (4) local minima of the energy functional can 
be defined only in an infinitesimally small region as a function of the polarization, the energy 
functional being perfectly fiat otherwise. 
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In order to address problems 1 and 2, following Nenciu, we restrict ourselves to periodic- 
polarized-insulating states, of which the lowest in energy is a long-lived resonance of the 
unrestricted system. Keeping this restriction in mind, we show that the position operator 
can be decomposed, in the crystal momentum representation, into a non-periodic part and a 
periodic part. The latter can be introduced in an ansatz Hamiltonian acting on the periodic 
part of the Bloch functions, from which the Berry phase formulation of the polarization is 
recovered, solving problem 3 as well. 

We are aware that this line of thought does not yet justify rigorously the use of this 
Hamiltonian: a more careful derivation, in the spirit of the mathematical work of Nenciu, 
would be needed. However, this rather simple Hamiltonian allows to recover all the previ- 
ously known lowest-order expressions for the polarization and its derivatives, and to derive 
other low-order expressions as well as generic expressions to all orders, as we shall see in the 
coming sections. 

Problem 4 is solved by introducing a regularization procedure in reciprocal space, sim- 
ilar in spirit to the real-space cutoff radius introduced by NV. For the regularized energy 
functional, the local minima have a finite basin of attraction as a function of the polarization. 

III. PERTURBATION THEORY APPLIED TO A VARIATIONAL 
TOTAL-ENERGY FUNCTIONAL 

In view of the application of perturbation theory to Eq. |10|, we summarize now the 
variational formulation of DFPT, as presented in Ref. |3^. We consider the formalism at 
its non-self-consistent level, without including the Hartree and exchange-correlation terms 
of the perturbative expressions. 

One considers a perturbative expansion of a variational principle applied to the electronic 
total-energy functional. In terms of the small parameter A associated with the perturbation, 
the perturbation series reads 

0{X) = + \o^^) + a2o(2) + A=^0(3) + ... , 
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(52) 

A=0 



for a generic observable O. The system Hamiltonian is H = K + Vext, and the total-energy 
functional is 

N 

E=J2 {^a\{K + Vext)\Va) , (53) 

where K and Vext are the kinetic-energy and external-potential operators. The total-energy 
functional is to be minimized under the orthonormality constraints for the occupied wave- 
functions 

{iPa \<^I3) = Saf3. (54) 

Using the Lagrange-multiplier method, the functional 

N N 

= (V^a + ^e^'*) I V^q) - [( V^a l^^p) " ^a/?] (55) 

a=l a,/3=l 

is minimized with respect to the wavefunctions. The minimum condition, 5F/5ip*^ = 0, leads 
to the Euler-Lagrange equation 

N 
/3=1 

Eq. ^ represents a set of generalized eigenvalue equations which assume the form of 
the usual eigenvalue equations when the so-called diagonal gauge is chosen to fix the phase 
arbitrariness of the wavefunctionsS Here, we keep the generalized form, as needed for the 
choice of gauge to be used in our theory. An expression for the Lagrange-multiplier matrix 
is obtained by multiplying Eq. ^ by an occupied wavefunction, leading to 

A/3a = {^p \H\ Lfa) . (57) 



We consider now the perturbation expansion of Eqs. p^-|57|. The orthonormalization 
condition becomes 

t{^i^Vr^) = o for ^>i- (58) 

i=o 
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The expansion of Eq. EB gives the generahzed Sternheimer equation 



j=0 



(i-J) 



N 



j=0 p=i 



(59) 



where if'-*-' = T^*^ + vl^J.f is the ith-ordei term in the expansion of the Hamiltonian. 
The expansion of the Lagrange-multipher matrix is given by 



j=0 k=0 



(k) 



(60) 



Finally, a generic term in the perturbative expansion of the total-energy functional in 
Eq. E^ is written 



N j i j 

a=l 1=0 k=Ol'=0 
N j j 

- E E E T.si^-i-k-nAfl{d^yp) 

a,f3=l 1=0 k=0 l'=0 



(61) 



where i = 2j or i = 2j + 1. We remark that only wavefunctions derivatives up to order A-^ 
appear in the ith-order term of the energy, as a result of the 2n + 1-theorem. Moreover, 
a minimum principle holds for E^"^^^ with respect to the jth-order variations of the wave 
functions, i.e., 5E^^^'^ /d^'i^ = 0. 



A particularly useful result derived in Ref. |3^ is a set of non-variational expressions for 
the second-order derivative of the energy. In the present work, the Hamiltonian is of first- 
order in the perturbation {v^^^^ = for i > 2), in which case the non-variational expressions 
are given by 
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The zeroth-order wave functions are chosen to obey the unperturbed eigenvalue equation 



by 



(0) 



AO) 



if^^^ ). From Eq. 16^, the zeroth-order Lagrange- multiplier matrix is given 



A (0) - 



(63) 
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In the present work, we use the so-called "parallel-transport" gauge, as discussed in 



Ref. |3^. In this gauge, the following condition is imposed on the derivatives of the wave 
functions 



(0) 



(64) 



which allows us to rewrite the expansion of the orthonormalization condition as 



(0) 



i-l 



(65) 



for 1 = 1 



IV. PERTURBATION THEORY APPLIED TO THE CONTINUOUS FORM 

A. Perturbation expansion and proof of gauge invariance 

Following the discussion in Sec. ||, we can develop a perturbation expansion for the 
electric-field problem. In this section, we discuss the continuous form of the theory. The 
cell-periodic Hamiltonian, including the perturbation term, is given in Eq. |13[ We apply the 
machinery of the variational DFPT to this Hamiltonian, by postulating that the expression 
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is to be minimized with respect to the {unk}, under the constraints 



(67) 



A local minimum will exist for the functional in Eq. ^ provided that a discretization of 
the fc-space integrals is performed. The continuum formulation which is considered in this 
section is valid only at infinitesimal fields. 
Applying Eq. we write 
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The unconstrained minimization of this functional is obtained by setting SF [{unk}] /^Unk 
0, leading to the corresponding Euler-Lagrange equation 



[H^' +ieS—\ \Unk) - J2 ^rnn{k) \Umk) = . 



(69) 



Next, we consider separately the perturbation expansions of Eqs. ^ and In both 
cases, we will demonstrate explicitly that the general expansion term transforms properly 
under a general unitary transformation of the occupied orbitals. 

1. Lagrange multipliers and orthonormalization constraints 



In the present case, Eq. |6^ for the orthonormalization constraints reads 
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, i = l 



giving the occupied-subspace projection of in terms of the lower-order solutions for the 
periodic functions. 

Since nj^^ = for alH > 2, and Hj^^ = ie^, the expansion of the Lagrange multipliers 



becomes 
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In the following development, we will make explicit use of the expressions for A^^(/c), 
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The second-order term reads 
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2. Energy 



The perturbation expansion for the energy is obtained from Eq. We analyze even 
and odd terms separately. For the even-order terms we write 
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while the odd terms are given by 
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An important aspect concerns the invariance of these expressions with respect to the 
choice of phases of the Bloch orbitals. More generally, we must consider unitary transfor- 
mations that keep the subspace of occupied states invariant. We show in Appendix ^that 
Eqs. |7^ and |7^ can be rewritten in such a way as to display the required gauge-invariance 
property explicitly. The lower-order derivatives are usually of more practical interest, and for 
that reason the invariant form of the energy terms up to fourth order are written explicitly 
here, along with the general expansion term. 

The second-order energy derivative is obtained by setting z = 1 in Eq. |73[ After some 
manipulation, this quantity can be written in the following form: 
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As expected, our formula for E^"^^ is identical to the linear-response expression.Ey 
The non-variational expression for E^'^^ (see Eq. ^) is given by 



2tt_ TV 

71 Jo \ ^ 



n=l 



4°' 



r(0) 
'nk 



U. 



27 



a 

2^ 



^22L N / ft N 

71=1 \ m=l 



(0)\ 
nk I 



U 



(0) 



^ ^ I (o)\ / (0) 



U 



(77) 



The fourth-order energy term is written 
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It is worth pointing out the simphcity of the expression for E'^'^^ ^ which mirrors that of 
almost exactly. 

The general even-order energy term for z > 2 is written 
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In the above and the following expressions, we use the notation [ie^ \ umk) (^mfc|) to 
indicate that d/dk acts only on the quantities embraced in parenthesis. In order to demon- 
strate that the energy derivatives fulfill the gauge- invariance requirement, we consider a 
general gauge transformationii among the occupied states at each k-point: 

\Unk) = E^'»"(^) l^"'^) ' (SO) 
m 

where U is an unitary transformation, i.e., UW = 1. It follows immediately that 
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Since ^ acts only on gauge-invariant quantities, Eqs. |7B| - |75| are themselves gauge- 
invariant. The same argument holds for the odd-order derivatives we derive below. 
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Turning now to the odd-order derivatives of the energy, we set i = 1 in Eq. |75| to write 
the third-order term as 

N / g N 
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This expression for E^^^ is identical to the one previously derived by Dal Corso and MauriS 
The general odd-order term for i > 1 is written 
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In appendix ^ we demonstrate how Eqs. [76| - [79| , p2| , and ^ are obtained from Eqs. ^ and 

m 



3. Sternheimer equation 



The projection of the wave functions on the subspace of occupied unperturbed states 
is given by Eq. The projection onto the subspace of unoccupied states is given by the 
projection of the Sternheimer equation in that subspace. The perturbation series for the 
Sternheimer equation can be obtained either by expanding Eq. ^ or more directly from 



Eq. |7| above, by setting /5 
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where Pck is the projector onto the subspace of unoccupied unperturbed states. This equa- 
tion can be solved for m^*^, once the lower-order derivatives of Unk and Am„(fc) have been 
obtained. 

Using the invariant form for the even terms of the energy, Eqs. [76| - |79| , we set 
^^(2«)^^^*W = to obtain the explicitly invariant form of the Sternheimer equation. For 
the i = 1 and i = 2 terms we obtain 
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For the complete specification of and u'^, the projections onto the unperturbed 
occupied subspace are obtained from Eq. [70|: 
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where P„fc is the projection operator for the occupied states. 

The higher-order terms for the Sternheimer equation are given by 
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The valence-band component of y is given by 
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As in the case of the energy terms, in Eqs. and the derivative d/dk acts only on 

gauge-invariant quantities. This completes our development of the perturbation expansion, 
and the proof of gauge invariance of the continuous formulation. 



B. Discretized form of lower-order expressions 

We examine now the discretized form of the lower-order terms for the energy and the 
Sternheimer equation. In practical calculations, it is mandatory to use a discrete set of 
/c-points to evaluate the Brillouin-zone integrals. However, when the focus is on E^'^\ the 
discretization of the ^ operation can be avoided, as the projection on the conduction bands 
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of the derivative of the wavefunctions versus k can be computed from a Sternheimer equa- 
tioni^l. This has the disadvantage to add a significant coding and computational step in the 
whole procedure. 

We choose the following symmetric finite-difference expansion for the derivatives with 
respect to k: 



dk l"-'^^^""^' ^ 2AA; 



(91) 



where Ak = kj^i — kj = {2TT/aNk). Clearly, this expression retains the gauge invariance of 
the continuous form. Next, Eq. ^is used in the derivation of explicit discretized expressions 
for the energy derivatives up to the fourth-order, and for the Sternheimer equation up to 
the second order. 

1. Energy 



From Eqs. |76| and pl| we obtain the discretized formula for E^'^'^: 
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The non-variational expression for E^'^'^ is written 
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The discretized versions of Eqs. ^ and ^ are 
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^. Sternheimer equation 



The discretized expressions for the i = 1 and i = 2 terms of the Sternheimer equation 
are given by 
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The required gauge invariance of the expansion terms is preserved in the discretized ex- 
pressions. Note that the solutions at a given k-point are now coupled to the first-neighbor 
k-points in the reciprocal-space grid. 



V. PERTURBATION THEORY APPLIED TO THE DISCRETIZED 

POLARIZATION 

A. General perturbation expansion 

Let us consider now the perturbation treatment of the problem on the basis of the energy 
functional given in Eq. |16|, where the polarization is written in a discretized form. This 
formulation can be viewed as the reciprocal-space analog of the NV@ real-space functional. 
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In this approach the gauge-invariance of the energy is guaranteed by the fact that Eq. |^ is 
itself gauge invariant .I'Ull 

We seek a minimum for Eq. |1^ with respect to the occupied orbitals under the 
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6mn- Lagrange multiphers are introduced to write the uncon- 
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In Appendix O, we prove the following result: 
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which can be used to derive the Euler-Lagrange equation from Eq. as follows: 
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Below, we consider the perturbation expansions of the Lagrange multipliers, the energy, 
and the Sternheimer equation. The expansion of the orthonormalization condition was 
already developed in the previous section, and remains unaltered in the present case. 

1. Lagrange multipliers 



We multiply Eq. |10(J| on the left by to write the Lagrange multipliers 
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It can be readily seen that the terms involving the overlap matrix cancel out, since 
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The perturbation expansion of Eq. |101| takes the simple form 
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2. Energy 



In order to write the perturbation expansion of Eq. we need the expansion of the 
polarization on the basis of the 2n+l theorem. The variation-perturbation framework allows 
us to focus on the part of E^^^^ or E^^j^^"* that comes from variation of the wavefunctions up 
to order i only. For these quantities, we introduce the notation E^^{^^ or E^^i^^'^\ The even 



terms are written 
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By the same token, for the odd terms we obtain 
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In these expressions, the 2n + 1-theorem implies that only the contributions of order 
< i from the perturbed wave functions will appear, when we consider the contribution of 
the polarization term to the total-energy derivatives. More explicit formulas for computing 
these polarization derivatives will be given below. 

With these results, we can expand Eq. From Eqs. ^ and |104| we obtain the even 
terms 
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while from Eqs. ^ and lUS, the odd terms are written 
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3. Sternheimer equation 



The perturbation expansion of Eq. |100| yields the Sternheimer equation 
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The conduction-band projector P^fe in Eqs- ^'^'^ |108| appears due the fact that the 
Lagrange multiphers in these two equations may be different, due to /c-gauge freedom. Thus, 
only conduction-band contributions can be identified while comparing the two formulations. 

By comparing Eqs. ^ and |108| , we see that in the present formulation the term 
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The theory of finite-difference approximations to derivatives could now be applied to 
Eq. |109| , as a function of Ak. Note that this expression is invariant under Ak —Ak, as it 
induces simultaneous exchange of fc^+i and Thus 

D{Ak) = D{0) + O{Akf. (110) 

One can now define 
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giving a higher-order approximation of -0(0) as 
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This improved expression also derives from a total energy functional. Instead of Eq. |T( 
one must start from 
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Despite its interest, we will not explore this topics further in the present paper. 
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B. Lower-order expressions 



In Appendix 0, we derive the Taylor expansion of Eq. |^ for the polarization, which 
allows us to obtain exphcit expressions for Eqs. |104| and |105[ Here, we look at the lower- 
order expressions for the energy and the Sternheimer equation. 



1. Energy 



From Appendix 0, the second-order polarization term is given by 
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where Q{kj, kj^i), obeying 
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is the inverse of the zeroth-order overlap matrix S^^{kj, fcj+i) 
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The second-order expression for the energy is then given by the i = 1 term in Eq. |10(j| 
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The third-order derivative of the energy is given by by the i = 1 in Eq. |107| . The first- 
order contribution to the Lagrange multipliers vanishes, since from Eqs. |7D| and |1U3| we have 
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which leads to E^'^'> = E^^i . From the results in Appendix 
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For the fourth- and higher-order energy derivatives, the expansion yields very involved 
expressions. We end this section by considering the fourth-order term for the energy in a 
more compact notation [we drop the {kj,kj^i) matrix arguments]: 
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In this expression, we write explicitly only the terms containing u^j^^, which will determine 
the second-order term of the Sternheimer equation. The corresponding fourth-order energy 
is given by 



Nk N 



k k=l n=l 



(0) (0) 



(4,2) 
pol 



(119) 



2. Sternheimer equation 
From Eq. p.08|, the first-order term for the Sternheimer equation reads 
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while the second-order derivative is written 
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These two expressions can also be consistently obtained by taking the conduction-band 
projection of the corresponding terms in the expansion of 6 E'^'^'^'' / 6u*^^. = 0. 

VI. NUMERICAL TESTS 

In this section, we illustrate the present theory by applying it the two-band ID-TB 
Hamiltonian introduced in Sec. |l|. For this model, exact analytical expressions can be writ- 
ten for the continuous formulation. Our purpose in this simple application is to demonstrate 
consistency between the continuous and the discretized versions of the theory, and also to 
make a preliminary assessment of the convergence of the energy derivatives obtained in the 
two discretized formulations, with respect to fc-point sampling in the Brillouin zone. Ap- 
pendix contains the detailed derivations of the results presented here. As such, it presents 
a step-by-step example of the use of the formalism developed in the present paper. 

A. Response to a homogeneous electric field 

In the continuous formulation, the first-order change of the valence state is obtained 



from the corresponding term in the Sternheimer equation. By setting i = 1 in Eq. |8J, this 
is given by 
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The solution to this equation is given in Appendix 0. 
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Having solved Eq. |T22|, E^^) 

can then be obtained from the simplest non-variational 

expression Eq. ^ as 
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where Ask = slk ~ ^vk^ dQk/dk is obtained from Eqs. ^ and ^ (see Eq. [Cil) . 

The second-order change of the valence state is obtained from the i = 2 term in Eq. |^ 
and from Eq. In our ID model, the Sternheimer equation for the conduction-band 
projection reads 
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Using the solution for m^^^ \ given in Appendix y, we can now arrive at the expression 
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for the fourth-order energy. Applying Eq. ITS], we write 
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(Expressions for E^'^^ and E^^^ in terms of elliptic integrals of the second kind are given in 
Appendix y.) 

Turning now to the discretized expressions, using the non-variational Eq. |7^ we compute 
the DAPE expression for the second-order energy: 
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The PEAD second-order energy expression takes a different form: 
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' i^^Asj \ 2Ak 

With the results in Appendix |C|, we have all the ingredients to write analytical ex- 
pressions for the discretized versions of E(^), but these are quite cumbersome and will not 
be reproduced here. The numerical results obtained using these expressions are discussed 
below. 
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B. Numerical results 



In order to test the consistency between the three formulations, we checked that by 
sufficiently increasing the number of /c-points used in the evaluation of the discretized energy 
derivatives, we obtain an agreement between continuous and discretized expressions within 
stringent degrees of accuracy. For example, for t = 1, 80 fc-points in the full Brillouin zone 
are needed for the three expressions for (Eqs. P|, and Wm to agree within ~ 1%, 



while 240 /c-points are needed to get the same level of agreement for the E^^^ expressions. 

By decreasing the number of /c-points, thus worsening the level of accuracy of the dis- 
cretized expressions, the differences between them become more apparent. In Fig. ^ we show 



the quantity -E^Jcr. ~ -^LLt I ^ exact giving the percentual error in the evaluation of E'^'^^ for 
the two discretized formulations, using 20 /c-points, with the hoping parameter varying over 



the [0,1] interval. In Fig. ^ the corresponding quantity for E''^\ E^-l^^ — E^^J^^^ / ElJ^^^f., is 
shown for a sampling of 80 /c-points. 

It is clear from Figs. ^ and ^ that the energy derivatives obtained from the PEAD 
converge faster with respect to the number of /c-points, at least for the present ID model. 
However, it must be borne in mind that this formulation involves the calculation of the 
inverse of the zero-field overlap matrix, and in practical applications the additional cost 
associated with this operation could offset the gain in fc-point convergence, specially when 
the two formulations are applied to systems with large numbers of atoms in the unit cell. 
This point remains to be further addressed when the theory is applied in the context of 
realistic tight-binding and ah initio calculations. 

We also computed the norm of the ffist- and second-order wavefunction derivatives 

1/2 

as a function of fc, for the continuous solutions and the two discretized 



forms, with a value of t = 1 for the hoping parameter and samplings of 20 and 80 /c-points. 
As expected from the above results for E'*^^^ and E'^^\ we observe that the wavefunctions in 
the PEAD are better approximations to the exact ones from the continuous formulation. 
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VII. SUMMARY 



The goal of this work was to obtain second- and higher-order derivatives of the total 
energy of periodic insulators, with respect to an applied homogeneous electric field. Related 
physical properties are the linear and non-linear dielectric susceptibilities (connected to 
linear and non-linear optical constants). 

Although a variation-perturbation framework had been formulated earlier for the compu- 
tation of derivatives of the total energy with respect to many different perturbations, several 
formal and technical difficulties must be addressed when considering the specific case of a 
homogeneous-electric-field perturbation. 

At the level of the electric-field-dependent energy functional, we proposed a basic expres- 
sion, Eq. that we argue to be valid in the space of states possessing a periodic density. It 
is directly linked to the modern theory of polarization, proposed by King-Smith and Vander- 
bilt nearly a decade ago, and allows to recover easily the Berry phase polarization formula, 
central to this theory. 

Unfortunately, when the polarization is varied, this energy functional leads to local min- 
ima with a basin of attraction of infinitesimal extent. A regularization procedure, based on 
a fc-point discretization of the reciprocal-space integrals, must be used in order to lead to a 
finite size basin. This is the reciprocal-space analog of the real-space cut-off introduced by 
Nunes and Vanderbilt in their treatment of polarized Wannier functions. 

Having thus defined a suitable energy functional, that depends on the applied homo- 
geneous electric field, we were allowed to proceed with the application of the variation- 
perturbation machinery. Interestingly, the derivation of the canonical formulas at all orders 
of perturbation can be done either on the basis of the energy functional already regularized, 
or on the basis of the unregularized one, followed by regularization at each order. The for- 
mulas derived in the two cases differ from each other. Working with the regularized energy 
functional gives more cumbersome expressions, however perfectly consistent with an energy 
functional, while the a posteriori application of the regularization at each order is not consis- 
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tent with the regularized energy functional. The two procedures will tend to the same limit 
when the discretization is refined further and further. The expression for the third-order 
derivative of the total energy, previously proposed by Mauri and Dal Corso, is recovered, as 
an instance of the " a posteriori" application of the regularization technique. 

We applied this formalism to a model one-dimensional two-band Hamiltonian, showing 
explicitly the pathology of a non-regularized energy functional and its cure, as well as the 
differences related to the order in which the perturbation expansion and the discretization 
procedure are applied. The two discretized formulations are shown to agree in the continuum 
limit, although the "perturbation expansion after discretization" (PEAD) formulation seems 
to be closer to the exact answer than the "discretization after perturbation expansion" 
(DAPE) formulation, for an equivalent grid. 
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APPENDIX A: 



In this appendix, we describe the algebraic manipulations needed to transform Eqs |7^, 
|75| , and into their gauge-invariant forms presented in Sees. |1V A 2| and [iV A 3[ We will 
use the following result: 
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Furthermore, from Eq. we have 
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With the help of Eq. ^ and 6, 
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it is straightforward to show that 
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We recall that the notation '\dk \ u) (m|)" indicates that dk acts only on the expression in 
parenthesis. 
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Now, for reference we repeat Eqs. ^ and |75| here: 
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Let us examine the i = 1 and i = 2 terms in Eq. |A4| , and the i = 1 term in Eq. We 
use Eq. ^ to write 
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Eq. ^ and the orthonormahty relation 
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two terms in this equation in the explicit gauge-invariant form of Eq. as follows: 
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where we use the notation \dkUmk) = ^ \umk)- The second term in Eq. ^ is obtained 
simply as the hermitian conjugate of this latter equation. 
For the third and fourth-order terms, we have 
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We now apply Eq. to recombine terms in these two expressions as follows: 
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where j = 1 applies to E^^^ and j = 2 to i?^^-*. This leads to the gauge-invariant expressions 
for these quantities, Eqs. ^ and |7|, respectively. 

For the general energy derivative, besides the terms corresponding to those above for 
E^^^ and E^^\ we must also consider terms of the form 
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For the second term on the right, we use Eqs. ^1| and |A2| to write 
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These results demonstrate the gauge invariance of the general expansion terms for the 
energy. The proper gauge-transformation properties of the Sternheimer equation can also 
be proven explicitly along the same lines, but it follows more simply from the invariance of 
the even-order terms for the energy. 



APPENDIX B: 



In order to develop the PEAD formulation, we examine the following expression appear- 
ing in Eqs. || and |l^: 

In det [Snm{,kj, kj+i){8)] = In det Sj^}^{kj, kj+i) + SSj^^{kj, kj+i) + ... 

= In det^S^^{kj,kj+i^ + dS—lndet[Snm{kj,kj+i){S)] , (Bl) 

where the perturbation expansion of Snm{kj, kj^i){S) is defined according to Eq. |52|. With 
the help of the "magic" formula 
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we can rewrite Eq. ^T] (to simplify the notation, in the following equations we drop the 
/c-point arguments) as follows: 
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To proceed further, we let AS{S) = S{S) - S^^^ = SS^^\S) + S^S^^\S) + S^S^^^{S) + 
and the inverse overlap matrix can be written as 

where, following the notation introduced in Sec. [VB| , we use Q = [S''^^]^^ for the inverse of 
the zeroth-order overlap matrix. 

To obtain the lower-order terms in the expansion of Eq. PT| , we write the terms up to 
third-order explicitly: 

S-\£) = Q- SQS^^^Q + [-QS^^^Q + QS^^^QS'^^^Q 

- Q5(^)Q5(^)Q5(^)q] + 0{8^) . (B5) 
Combining Eqs. il, il and dS{S)/dS = S^^^ + 2SS^^^ + 38"^ S^^^ + 48^ S^^^ + 0{S^), we 



arrive at 
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Finally, we note that the magic formula may also be used to derive the Euler-Lagrange 
equation given in Eq. |100| , as follows: 
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APPENDIX C: 



In this appendix, we provide all the steps for the derivation of the energy derivatives for 
the ID-TB model, given in Sec. [VI| . The expressions for the discretized versions of are 
not written explicitly. 

Eq. ^ for the unperturbed total energy is a complete elliptic integral of the second 
kind§, which is given in its general form by 

JJ dy {l + A' cos' yy , (CI) 

where n is a positive or negative integer, and A = 4t. Several such integrals will be encoun- 
tered in the course of our derivation. 



1. Continuous formulation 



We begin by developing some useful preliminary results. For our ID-TB model, the 
derivatives with respect to k of the unperturbed valence and conduction states (Eq. ^) are 
written analytically as 
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where Aa^ = a^k - OLvk- 
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where Asu = e 
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Using Eq. the second-order energy term (Eq. |123D is written 
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The second-order change of the valence state is obtained from Eqs. ^ and The 
first-order Lagrange multipher in Eq. |86| is obtained from Eqs. p8|, |72|, and 
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Simple algebraic manipulations involving Eqs. |C3| , and [061 , combined with Eqs. ^ and 
|88|, yield the second-order wave-function derivatives 
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These results lead to Eq. |125| for E'^^\ In terms of elliptic integrals, this quantity is 
written 
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2. DAPE formulation 



We now apply the DAPE expressions obtained in Sec. |IVB| to our ID-TB model. The 
first-order valence state is given by 
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Using Eq. BHl this expression becomes 



(1) 

vki 



2AkAe^ 



sin (6j+i — 0j-i) cos (20j — Qj+i — ©j-i) e 



U 



(0) 



(CIO) 



51 



where we use a simplified notation 0^^ — > Qj (likewise for Aa and Ae). From Eqs. piO| and 
93|, we obtain Eq. |I2^ . 



The second-order wavefunction is obtained from Eqs. ^ and |9^. The Sternheimer equa- 
tion is written 
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By combining Eqs. p8|, pSl and piO|, we arrive at 
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Eqs. CIO and C12, combined with Eq. 95, lead to an analytic expression for E^^\ 



3. PEAD formulation 



Here, we apply the PEAD expressions from Sec. |V| to the model system. The zero-field 
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Note that for this model, S" is a 1x1 matrix, and hence the inverse overlap matrix is simply 
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From Eq. |120|, the first-order Sternheimer equation reads 
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Using Eq. 38, we obtain 
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Plugging the above results for Q, 5*°, 5"^^^ and m^,],-* in Eqs. |8] and |121| we get 
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As in the DAPE case, from these expressions for u'^l, and u~^^. follows the PEAD analytic 
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FIGURES 

FIG. 1. The phase in the trial wavefunctions for the ID model. 7^ changes by 2tt over a 
small interval A/c centered at an arbitrary A;-vector in the Brillouin zone. The position of the jump, 
{k), is indicated in the figure by the vertical dotted line. 



FIG. 2. The function cos7jfc, which differs from 1 over an interval AA; in the Brillouin zone. 



FIG. 3. The percentual error 
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exact 



in the evaluation of the second-order 



change in energy for ID model using the DAPE and the PEAD formulations, with a 20 A;-point 
sampling of the Brillouin zone. The hoping parameter varies over the [0,1] interval. 



FIG. 4. The percentual error 
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exact 



in the evaluation of the second-order 



change in energy for ID model using the DAPE and the PEAD formulations, with an 80 A;-point 
sampling of the Brillouin zone. The hoping parameter varies over the [0,1] interval. 
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